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ABSTRACT 

Computations of cosmic-ray transport based upon 
finite-diference methods are afflicted by instabilities, 
inaccuracies, and artifacts. To avoid these problems, we have 
developed a Monte Carlo formulation which is closely related not 
only to the finite-difference formulation, but also to the 
underlying physics of transport phenomena. Implementations of 
this approach are currently running on the Massively Parallel 
Processor at Goddard Space Flight Center, whose enormous 
computing power overcomes the poor statistical accuracy that 
usually limits the use of stochastic methods. These simulations 
have progressed to a stage where they provide a useful and 
realistic picture of solar energetic particle propagation in 
interplanetary space. 

Keywords: Cosmic-rays, Particle Transport, Interplanetary 
Medium. 


INTRODUCTION 


Over the past three years the computations have evolved through 
simulations of rectilinear transport along a constant guiding field, 
which were presented at the First Symposium on Massively 
Parallel Scientific Computation (Ref. 5), and of focused 
transport with constant focusing length and mean free path (Ref. 
6), to the present simulations, which allow arbitrary variations of 
the focusing length and mean free path. In the interplanetary 
context, this formulation includes all important effects except 
those of convective motion of the background medium, which 
significantly affect the slow variations of cosmic-ray modulation, 
but which play a minor role during the rapid evolution of solar 
particle events. It includes the two essential aspects of charged 
particle transport. These are a strong inhibition of transport 
perpendicular to the guiding field and a strong anisotropy of the 
pitch-angle scattering by random fields. 

TRANSPORT EQUATIONS 

Under the circumstances outlined above, particle transport is 
described by 


The diffusion idealization, which has been almost universally 
invoked in discussions of cosmic-ray transport is easy to treat 
analytically. However, many observed phenomena give clear 
evidence for non-diffusive effects. One example is the so-called 
"scatter free" propagation of kilovolt solar electrons (Ref. 1), 
which is inconsistent with diffusion, but which can readily be 
interpreted in terms of a coherent mode of propagation. This 
mode is novel, but it is just a manifestatation in a dynamic 
situation of non-diffusive effects similar to those considered in 
the steady-state by classical transport theory (Ref. 2). Although 
these effects have been described analytically in References 3 and 
4, the theory is very complicated. Consequently, there is a need 
for reliable numerical computations which bypass these 
complexities and yield concrete results suitable for comparison 
with observations. This paper describes computations of charged 
particle transport along a large-scale guiding magnetic field B 
whose spatial variations are characterized by the tocusing length 
L, which is defined by 
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in which h is the number of particles per unit distance parallel to 
the guiding field, p is the cosine of the pitch-angle. The 
parameter s = Vt, where V is particle velocity, plays the role of a 
temporal variable. The Fokker-Planck coefficient of pitch-angle 
scattering is given by 
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where q is an index that measures the anisotropy of scattering 
(Ref. 7). The function G that appears in Equation 2 is defined in 

terms of by 


where z is distance parallel to the guiding field. Particle trajec- 
tories are scattered by small-scale random magnetic fields whose 

effect is described by the mean free path K. Note that the 
magnetic fields are visualized as static and that there is no 
interaction among particles in an extremely tenuous distribution 
of charged particles. This situation differs from those considered 
by plasma physics, but it is closely analogous to those treated by 
classical transport theory. 
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In the discrete formulation, the continuous variables are replaced 
by a three-dimensional grid whose spacings are A z, Ajli, and As, 
and the derivatives appearing in equation (1) are replaced by their 
finite -difference analogs. The effect of these replacements is best 
described in terms of the particles flowing in and out of a cell, 

whose dimensions are Ap and Az, during a temporal increment 

As. These flows are illustrated in Figure 1, where the integer M 
refers to the pitch angle, and the integer Z refers to the distance z. 
Here, particles passing through the top and bottom of the box are 
those whose pitch angle changes because of scattering, and 
particles passing through the sides are those whose distance 
changes as a result of their motion parallel to the guiding field. 
The flows due to scattering can be described in terms of 
coefficients pM and qM, which characterize the flows out of box 

M toward larger and smaller values of p, respectively, and 
whose numerical values can be derived from Equation 2. The 
difference equations that correspond to Figure 1 can be solved by 
standard numerical methods, but these methods are difficult to 
implement and subject to subtle errors, which are discussed in 
Reference 8. 



Figure 1. The number of particles in each cell of the p - z plane 
changes, because scattering causes particles to flow in and out 
through the top and bottom boundaries, while motions in space 
cause them to flow through the vertical txmndaries. 


In the Monte Carlo formulation, the random history of a large 
number of particles is followed under the assumption that the 
coefficients p and q can be interpreted as probabilities in each 

temporal step that p will change by Ap toward more forward or 
more backward directions, respectively. Obviously, this defines 
1-p-q as the probability that a pitch angle will not change. After 
their pitch angles are updated in each step, the particles move a 
distance p As. 

In Figure 2, the transition probabilities are plotted as a function 


of p for strongly anisotropic scattering, q = 1.8, similar to that 
occuring in interplanetary space and for strong focusing, = 5. 
These probabilities exhibit the same gross features as the Fokker - 
Planck coefficient, but there is a significant assymetry such that p 
is consistently larger than q. This leads to the systematic drift 
toward forward directions that is expected as a result of focusing. 

Note that the probability of going backward through p = 0 is 
only 0.13% in this example. This leads to a coherent mode of 
transport that is very different from diffusion, for any particle 
that reaches the forward hemisphere has a very small chance of 
going backward again. 



Figure 2. Forward and backward transition probabilites plotted 
as a function of pitch-angle cosine. 


THE ALGORITHM 

To implement the Monte Carlo scheme outlined above, each 
particle was assigned an integer distance and an index M cor- 
responding to p that lies between 0 and 8. Because these 
parameters occupy only three bytes for each particle, there was 
plenty of storage for several parallel arrays of particles. 
Consequently, the results given below are based on 32 arrays 
which contained 32 * 16384 = 524288 particles. The fate of the 
particles was determined by a single parallel array of random 
integers ( ranging from -32767 to +32767) that was updated 
frequently. To implement changes in the pitch-angle cosine, 
forward and backward integers were assigned to each particle, 
according to its pitch angle and distance, in such a way that the 
probability of the current random number being larger than the 
integer is the corresponding forward or backward probability 
analogous to those poltted in Figure 2. Then the angular index 
was incremented for those particles whose current random 
integer was positive and greater than the forward integer, and 
decremented for those whose random integer was negative and 
less than the backward integer with its sign reversed. This 
approach satisfies the basic requirement that the probabilities of 
incrementing, decrementing and leaving unchanged the 
pitch-angle must add to unity. After the pitch-angles had been 
updated, each particle’s distance was incremented. When the 
desired number of temporal steps had been carried out, particles 
were binned according to distance and pitch-angle. 
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Figure 3. Density profile for rectilinear transport after a 
collimated injection at the arrow. 

SPATIAL AND TEMPORAL PROFILES 

To set the stage for the simulations of interplanetary transport that 
are the main subject of this paper, it is useful to consider tirst the 
case of rectilinear transport, in which the guiding field is 
constant, and focusing does not occur. Figure 3 presents results 
obtained from the MPP as plots of the total number of particles in 
each distance bin. This sum over pitch-angles is a measure of the 
isotropic particle density. Because the total number of particles 
was large, statistical errors are small and, consequently, are not 
shown explicitly. However, slight irregularities in some parts of 
the curves give an indication of their magnitude. 

This density profile describes a situation very shortly after the 
injection of a collimated beam of when the particles have had time 
to move a maximum distance of only one mean free path. The 
curve exhibits two features: an intense localized peak at the right 
and a broad wake spread on both sides of the arrow which 
indicates the point of injection. Qualitatively, the peak appears 
because particles become nearly uniformly distributed in the 
forward hemisphere, while very few particles penetrate to the 
backward hemisphere through the region of weak scattering at 

(ji = 0. This means that the particles in the forward hemisphere 
move with nearly the same average velocity parallel to the field, 
but statistical fluctuations in individual velocities give rise to a 
peak centered around the average displacement. Such features 
are designated as coherent pulses. They decay exponentially as 
particles slowly escape into the backward hemisphere and join 
the wake. 

In the spiral pattern of the interplanetary magnetic field, the 
focusing assymetry that appears in Figure 2 tends to keep 
particles in the forward hemisphere. Consequently, the coherent 
pulse is long lived and intense relative to the wake. Moreover, 
the assymetry becomes very pronounced near the sun, for the 

ratio X/L varies approximately as 1/r, where r is distance from the 
sun. 

To illustrate these effects, Figure 4 shows snapshot profiles of 
number of particles vs. r for four equally spaced times after 
injection into a model interplanetary field at r = 0. These profiles 


describe a moving pulse whose width increases as it moves out, 
which is the qualitative behavior expected. However, the exact 
evolution of this width and the deceleration of the peak, which is 
evident from a careful examination of the figure, are details 
which crucially affect observations, but which are not adequately 
described by the analytic theory of focused transport. 

In contrast with the rectilinear profile illustrated by Figure 2, the 
wake is virtually invisible in all of the four profiles. As was 
discussed above, this enhancement of the coherent mode is a 
consequence of focusing. 



Figure 4. Profiles of isotropic particle density vs. distance from 
the sun for four different times after injection. 



Observations from spacecraft generally take the form of profiles 
of intensity vs. time at an essentially fixed position in space. In 
Figure 5, results from MPP simulations, which are presented in 
this form, show what happens when the coherent pulse sweeps 
over an observer at a fixed location. The dotted curve, which 
refers to a mean free path equal to the distance of the observer 
from the sun, defines a peak followed by a very small wake. As 
the mean free path becomes smaller, this pulse becomes less 
prominent relative to the wake ( solid curve ), until it is 
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submerged in the onset ( dashed curve ). Examples of all three 
of these behaviors are seen among observed solar particle events, 
along with events whose extremely slow onset can be understood 
in terms of standard diffusion theory. Consequently, the MPP 
results make it clear that all types of prompt solar particle events, 
diffusive, abrupt onset and coherent, can be understood as 
morphological stages in a continuous sequence of profiles 
predicted for increasing mean free paths by the theory of charged 
particle transport. The same basic physics explains all types. 

ANGULAR DISTRIBUTIONS 

The discussion above has focused on isotropic intensities in 
which particles are counted regardless of their directions. In fact 
both predicted and observed intensities are often highly 
anisotropic, and valuable information can be derived from a 
careful analysis of these anisotorpies. To illustrate this point, 

Figure 6 presents numbers of particles binned according to p and 
summed over a coherent pulse. Evidently, the intensity in the 
forward hemisphere p > 0 is overwhelming larger than that in the 
backward hemisphere. The solid curve which gives the result of 
a MPP simulation is in very good agreement with the dashed 
curve, which gives the function exp ( G ) which is expected from 
theory. 
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Figure 5. Angular distributions for particles summed over the 
coherent pulse. 
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CONCLUSIONS 

Results obtained on the MPP with the aid of Monte Carlo 
methods are equivalent in every detail to those based upon careful 
use of more traditional methods, but they are less subject to error 
and are closer to the physics. These characteristics offer 
tremendous advantages in the investigation of exotic transport 
regimes for which no theoretical description is available. In 
particular, the formulation of problems in which particles gain or 
lose energy leads to prohibitively large conventional com- 
putations, but their Monte Carlo versions are not significantly 
more complicated than the one described here. We intend to 
exploit these advantages in the investigation of two such 
problems: Adiabatic deceleration of cosmic-rays due to 

expansion of the solar-wind, and the loss of energy by electrons 
in radio sources due to synchrotron radiation. On a more 
immediate time scale, the current simulations, which embody an 
accurate and useful description of solar particle events, will be 
applied to the interpretation of observations. 
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